function compute_DWL_bootstrap(mat_type,range_extend,K1,K2,K3,trim_perc_q,y_min,y_max,constr,RR,berkson_txt,factor_berkson,berkson_off,input_directory,output_directory,fileID)

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% COMPUTE BASELINE ESTIMATE
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Berkson_estimated=0;
unit_TT=0;

[data,config]=berk_ml.load_data(mat_type,range_extend,K1,K2,K3,trim_perc_q,y_min,y_max,berkson_txt,factor_berkson,berkson_off,Berkson_estimated,unit_TT,input_directory,output_directory);

estimates=berk_ml.estimate(data,config,constr);
q=data.q;

y_interest_univ_DWL=config.y_interest_univ_DWL;
[DWL_pertax,DWL_perinc]=get_DWL(estimates); 

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% SET SEED
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
rng(123456,'twister');

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% BOOTSTRAP LOOP
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
DWL_pertax_bootstrap={};
DWL_perinc_bootstrap={};

data_bootstrap={};
config_bootstrap={};

for rr=1:RR

	[data_bootstrap{rr},config_bootstrap{rr}]=berk_ml.load_data_bootstrap(mat_type,range_extend,K1,K2,K3,trim_perc_q,y_min,y_max,berkson_txt,factor_berkson,berkson_off,Berkson_estimated,unit_TT,input_directory,output_directory);

end

parfor rr=1:RR
	
	estimates_bootstrap=berk_ml.estimate(data_bootstrap{rr},config_bootstrap{rr},constr);
	[DWL_pertax_bootstrap{rr},DWL_perinc_bootstrap{rr}]=get_DWL(estimates_bootstrap);

end

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% EVALUATE
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
N=size(q,1);
DWL_pertax_bootstrap_vec=NaN*ones(RR,1);
DWL_perinc_bootstrap_vec=NaN*ones(RR,1);

alpha=10;

	fprintf(fileID,['\n --- CI AT ' num2str(alpha) '%% LEVEL ---\n']);

	for tt=3:(-1):1
		for yy=1:3
			inc=exp(y_interest_univ_DWL(yy));
			DWL_pertax_baseline=DWL_pertax(tt,yy);
			DWL_perinc_baseline=DWL_perinc(tt,yy);
			for rr=1:RR
				DWL_pertax_bootstrap_vec(rr)=DWL_pertax_bootstrap{rr}(tt,yy);
				DWL_perinc_bootstrap_vec(rr)=DWL_perinc_bootstrap{rr}(tt,yy);
			end
			[CI_low_pertax, CI_high_pertax] = get_ci (DWL_pertax_baseline, DWL_pertax_bootstrap_vec, alpha, 1, N, N, 'equal');
			[CI_low_perinc, CI_high_perinc] = get_ci (DWL_perinc_baseline, DWL_perinc_bootstrap_vec, alpha, 1, N, N, 'equal');
			fprintf(fileID,'tt=%3.0f, yy=%3.0f, inc=%10.0f, DWL-tax=%10.3f [%9.3f; %9.3f], DWL-inc (x 10^4)=%10.2f [%9.2f; %9.2f]\n',...
				tt,yy,inc,DWL_pertax_baseline,CI_low_pertax,CI_high_pertax,10^4*DWL_perinc_baseline,10^4*CI_low_perinc,10^4*CI_high_perinc);
		end
	end

	%

